Quantum Monte Carlo with Coupled-Cluster wave functions 
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We introduce a novel many body method which combines two powerful many body techniques, 
viz., quantum Monte Carlo and coupled cluster theory. Coupled cluster wave functions are intro- 
duced as importance functions in a Monte Carlo method designed for the configuration interaction 
framework to provide rigorous upper bounds to the ground state energy. We benchmark our method 
on the homogeneous electron gas in momentum space. The importance function used is the coupled 
cluster doubles wave function. We show that the computational resources required in our method 
scale polynomially with system size. Our energy upper bounds are in very good agreement with 
previous calculations of similar accuracy, and they can be systematically improved by including 
higher order excitations in the coupled cluster wave function. 



Introduction. — Quantum Monte Carlo (QMC) meth- 
ods have now become standard tools for computations 
in a wide variety of strongly correlated systems [l| rang- 
ing from quantum chemistry [2[ to condensed matter [3| , 
and nuclear physics [4|]. The most celebrated capability 
of such methods is to provide very accurate estimates 
of ground state and thermodynamic properties of many- 
body systems, while enjoying a very favorable scaling of 
computational time with system size. 

In fermionic systems, QMC methods are plagued by 
the 'sign problem'. In principle, the sign problem can 
be circumvented in the 'fixed-node' (FN) approximation 
with the help of importance sampling with a trial ground- 
state wave function [5[ . The spectacular success of the 
FN QMC is largely because of the development of high 
quality trial ground-state wave functions. However, most 
of these wave functions have forms which are convenient 
for calculations in coordinate space (r-space). 

The lack of accurate and computationally efficient trial 
ground-state wave functions has, thus far, precluded a 
wide exploration of these algorithms within the config- 
uration interaction (CI) scheme (see, however, Ref. m). 
Generally, QMC methods within the CI scheme tend 
to rely on auxiliary fields introduced via the Hubbard- 
Stratonovich transformation [7]. Much recent interest, 
however, was sparked by the demonstration that even 
within the CI scheme it is possible to apply stochastic 
projection to systems much larger than what would be 
possible using conventional matrix diagonalization 8] . 

The development of an efficient ground state QMC al- 
gorithm within the CI scheme/momentum space would 
be of great interest in nuclear physics, where most of the 
modern interactions are written in non-local forms [9j; 
and in electronic structure calculations with non-local 
pseudopotentials [id ]. 

On the other hand, coupled cluster (CC) wave func- 
tions provide very accurate and size extensive approxima- 
tions for the ground state wave function for many phys- 
ical systems of interest Uj. The CC energies towards 



accurate QMC calculations in many-body systems [12J . 
In fact, CC calculations are considered to be the 'gold 
standard' in quantum chemistry [13j . More recently, CC 
theory has also reemerged as a method of choice in nu- 
clear structure calculations 
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In this Letter, we introduce a novel scheme to combine 
these two powerful many-body techniques. We propose 
to use CC wave functions as importance functions for a 
projection quantum Monte Carlo algorithm within the 
configuration interaction scheme, which we call config- 
uration interaction Monte Carlo (CIMC). In CIMC the 
ground state wave function of a CI Hamiltonian is fil- 
tered out by propagating the amplitudes of an initial 
arbitrary wavefunction via a random walk in the many 
body Hilbert space spanned by a basis of Slater deter- 
minants [15[. The CC wave functions are used to guide 
this random walk via importance sampling in order to 
circumvent the sign problem. Our method provides rig- 
orous upper bounds to the ground state energy whose 
tightness can be systematically improved by including 
higher order excitations in the CC wave function. 

We apply our method to the three dimensional homo- 
geneous electron gas (3DEG) in momentum space. The 
3DEG is described by a simple Hamiltonian; it neverthe- 
less encapsulates many of the difficulties associated with 
modern many-body theories. In particular, it has both 
the weakly and strongly correlated regimes which can 
be accessed via a single tunable density parameter, the 
Wigner-Seitz radius r s , thus providing an ideal system 
for benchmaking many body theories 16H21|. 

Monte Carlo in Fock space. — The CIMC projection 
algorithm is discussed in detail in Refill We describe 
it briefly for completeness. Consider a general second 
quantized fermionic CI Hamiltonian which includes only 
two-body interactions (atomic units will be used through- 
out) 
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where a] creates a particle in the single-particle (sp) state 
labeled by i (i is a collective label for all sp quantum 
numbers). The set S of sp states is assumed to be finite 
and of size Af s ■ The Vff are general two-body interaction 
matrix elements. 

For homogeneous systems we can use the plane wave 
states, with definite momentum and spin, as the sp basis 
set. The sp energies are £j = kf/2m, where k^ is the 
momentum of the i-state, and m is the fermion mass. 
We include in our sp basis all single-particle states i with 
kf < fc„ lax - In principle, cutoff independent results can 
be obtained by performing successive calculations with 
increasing fc max and then extrapolating to fc max — >• oo. 



For the 3DEG, the V°j b are given by 



Vij — (1 - < 5k a -k,.o)^k,+k j ,k a +k i; 
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The volume $7 of the simulation box (and hence the unit 
spacing in momentum) is determined by the density and 
the number of particles N in the simulation. We ignore 
the Madelung term because it does not affect the corre- 
lation energy. 

The many-body Hilbert space is spanned by the set of 
all .ZV-particle Slater determinants constructed from the 
sp orbitals i e S. We will denote these Slater determi- 
nants or 'configurations' with the vector notation, |n), 
where n = {n,} and rii = 0, 1 is the occupation number 
of the sp orbital i in | n) . 

The ground state wave function fy gs of H can be pro- 
jected out by repeated applications of the projection op- 
erator V = I-At(H-Et), on some initial wave function 
^o which has a non-zero overlap with ^ gs , 
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where the energy shift Et and the imaginary time step 
At are related by At < 2/(E max — Et), with _E max being 
the maximal eigenvalue of H [22( . 

In a Monte Carlo algorithm the wave function at any 
time step M, \&m, is represented as an ensemble of con- 
figurations. The one time-step projection (single appli- 
cation of V) is performed independently and stochasti- 
cally for each configuration, by interpreting the matrix 
elements (m|7'|n) as probabilites. 

In general (the 3DEG inclusive), V will have one 
or more negative off-diagonal matrix elements, which 
cannot be interpreted as probabilities. In principle, a 
stochastic evolution can still be carried out by evolv- 
ing 'signed' configurations, but this leads to an expo- 
nential decay in the signal to noise ratio with increasing 
time-step M . This is a manifestation of the sign prob- 
lem. Recently, it was shown that the sign problem can 
be somewhat moderated by including a configuration- 
annihilation step in the evolution [23j. Nevertheless, even 
the latter algorithm has exponential scaling with the sys- 
tem size, though with a reduced exponent. 



In CIMC, we circumvent the sign problem with the 
help of an importance function <I> G . We define a family 
of Hamiltonians, Jty, whose off-diagonal matrix elements 
are given by [2J, |25j , 
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while the diagonal matrix elements are given by, 

<n|H 7 |n) = (n|ff|n) + (l + 7 ) ]T s(m,n). (5) 

s(m,n)>0 

where, s(m, n) = <I> G (m)(m|_ff|n}/$ G (n). 

In addition, we define a new propagator V 1 as, 

<m|7> 7 |n) = 1 - Ar$ G (m)(m|H 7 - £ T |n)/$ G (n) . (6) 

The propagator "P 7 , by construction, is free from the 
sign problem for 7 > 0, and filters out the wave func- 
tion <I> G (n)\I/ 7 (n), where ^(n) is the ground state wave 
function of T~i 1 . 

The ground state energy £ 7 of H 1 is an strict upper 
bound for the ground state energy E gs of the true Hamil- 
tonian H for 7 > 0, and this upper bound is tighter than 



the variational upper bound ($ G |i/|$ G ) [la |2J, |25( . In 
addition, a linear extrapolation of £ 7 from any two values 
of 7 to 7 = — 1 also provides a rigorous upper bound on 
Egs [26( | . We found the extrapolation using 7 = 0, 1 to 
be best compromise between accuracy and statistical er- 
ror. Thus, the tightest upper bound for the ground state 
energy E is -Ecimc = 2£ 7=0 - £ 7 =i. 

The simple projection algorithm described above be- 
comes extremely inefficient for large -E max , i.e, for large 
N or jV" s - In practice, we use a much more efficient algo- 
rithm (free from any time-step error) which was proposed 
in Ref. El (see also, Ref. [H). 

Coupled cluster importance functions. — A good choice 
for the importance function $ G is crucial for the success 
of our method. We need <& G to be sufficiently flexible to 
include the dominant correlations in the system, and yet 
be quick to evaluate. In many strongly correlated sys- 
tems, CC wave functions fulfill the first criterion. Below, 
we describe a recursive algorithm which can be used to 
evaluate CC wave functions very efficiently 

The CC wave function can be written as 



|$cc) = e T |$ ) 
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where $0 is a model wave function and T is an op- 
erator which generates excitations on <I>o. We choose 
^0 = 'I'hf, i-c, the Hartree-Fock wave function. Due 
to momentum conservation, in a homogeneous system 
the simplest non-trivial approximation for T is T = T2 — 
^2ij ab tij a a a b a j a i' wnere s p states z and j (a and b) are 
occupied (unoccupied) in $hf- This constitutes the so- 
called coupled cluster doubles (CCD) approximation for 
the ground state wave function. 



Writing the wave function in terms of the excitations 



on top of $hf it can be shown that the only non-trivial 
components of the wave function are given by 
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for even to > 0. In the above equa- 

tion, *8bD(EK™) = *cc D (n) for |n) = 



,t 



■■,,, Pm O/n ...o/, m |$HF)) with pi < ... p m 

and hi < . . . < h m . The wave function has a vanishing 
component for odd to, and the component for to = is 
hxed by our choice of normalization to be 1. 

Using Eq. ([5]) we can easily calculate 3>ccd recur- 
sively. The computational effort for a single calculation 
has a combinatorial scaling with the average number of 
particle-hole excitations in $<3. However, it does not de- 
pend on N or J\f s ■ 

Although for this exploratory study we have used the 
CCD wave function, it is evident that similar recursive re- 
lations can be easily written for other CC type wave func- 
tions, e.g., CCSD or CCSDT, which will become neces- 
sary for applying our method to inhomogeneous systems 
or, systems with stronger correlations or many-body in- 
teractions. 

Results for the homogeneous electron gas. — In Table U 
we show the CCD energies, calculated using conventional 
CC theory, along with the corresponding Monte Carlo 
energies of an 3DEG system with N - - 14 and Af s = 342, 
for r s = 0.5, 1.0 and 2.0. We see that, for r s = 0.5 the 
CCD energy is very close to the corresponding Monte 
Carlo energy. But, for r s — 1.0 and 2.0 they are, in 
fact, lower than the corresponding Monte Carlo energies. 
Since, -Ecimc < (^ccdI-H^ccd) (see discussion above), 
this shows, once again, the non-variational nature of the 
energies obtained from conventional CC theory |27J |. 



Correlation energy (a.u.] 
CCD + CIMC CCD(l) ' 



CIMC 



0.5 

1.0 

2.0 



-0.572682 
-0.506701 
-0.417946 



-0.5729(3) 
-0.5021(3) 
-0.40317(2) 



-0.659641 
-0.657347 
-0.665071 



-0.5733(2) 
-0.5025(2) 
-0.4029(3) 



TABLE I. Correlation energies for AT = 14 and Af s = 342 
from conventional CC theory with the CCD and CCD(l) wave 
functions, along with the corresponding CIMC energies using 
each as importance functions. The numbers in parenthesis 
indicate statistical error in the last significant digit. 
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demanding ways of obtaining the i" b , while still preserv- 
ing the structure of the CC wave function. 

A simple option is to compute tfj by the second order 
M0ller-Plesset perturbation theory (MP2). If the CC 
equations are solved iteratively, then this is equivalent 
to stopping after the first iteration. In Table U we also 
compare the energies obtained from our CIMC projection 
using the MP2, denoted by CCD(l). 

The CCD(l) amplitudes produce a worse approxima- 
tion to the ground state wave function as compared to the 
full CCD amplitudes. Nevertheless, when used as impor- 
tance functions in our CIMC algorithm, the final estimate 
for the ground state energy for both cases are very close. 
The statistical errors are comparable for r s — 0.5 and 
1.0. For r s = 2.0 they are about an order of magnitude 
lower when the full CCD amplitudes are used. 

The above observation is extremely encouraging be- 
cause it means that, for our purposes, it may not nec- 
essary to solve the full CC equations to get reasonable 
tfj amplitudes. Of course, we do not expect the CCD(l) 
amplitudes to be satisfactory for more strongly corre- 
lated systems. Still, even in those cases one can, presum- 
ably, use computationally inexpensive approximations to 
the full CC equations. For the rest of this work, all the 
CIMC results have been computed using the CCD(l) am- 
plitudes. 

In Fig. [1] we show the CIMC ground state energy esti- 
mates for N — - 14 and r s = 0.5, 1.0, 2.0 and 3.0 for some 
of our large basis size calculations. In Refs. 120, |2lj, and 
28l it was suggested that for the 3DEG it might be pos- 
sible to extrapolate to the J\f s — > oo limit by exploiting 
a linear 1/Af s dependence of the correlation energy for 
large but finite Af s . Although, for r s — 0.5 and N = 14 
such a linear trend in the correlation energy is visible, 
for the other values of r s shown in the figure, no such 
trend is evident. The situation is similar for calculations 
we performed with N = 32 and 54. Thus, at least up to 
our largest basis size Af s = 2378, we cannot safely do an 
extrapolation to Af s —> oo with a reasonably low x ■ 

In Table Ull we show the ground state energy of different 
r s and N for the largest Af s (= 2378) calculated by us. For 
comparison, we also include the energy estimates from 



In principle, the CC amplitudes tf] can be obtained Refs. [IJ and |21 



from conventional CC theory. However, solving the CC 
equations is computationally expensive. Therefore, we 
investigated the possibilty of using less computationally 



The energies in Ref. [l9| are calculated using the r- 
space diffusion Monte Carlo method with an importance 
function that included backflow correlation on top of the 
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-0.55 
-0.56 
-0.57 
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FIG. 1. Correlation energies for N — 14 and r 3 — 0.5, 1.0, 2.0 and 3.0 as a function of the single particle basis size from CIMC 
with the CCD(l) importance function. The statistical errors are smaller than the size of the symbols. The lines are drawn as 
a guide to the eye. 







Correlation 


energy (a.u) 


r a 


N 


CIMC 


Other 


0.5 


14 

38 
54 


-0.5875(6) 

-1.809(4) 

-2.354(2) 


-0.5959(7) [21] 
-1.849(1) [2TJ 
-2.435(7) [21] 
-2.387(2) [19J 


1.0 


14 

38 
54 


-0.5114(5) 

-1.521(4) 

-2.053(4) 


-0.5316(4) [21] 
-1.590(1) [21] 
-2.124(3) [21] 
-2.125(2) [19] 


2.0 
3.0 


14 

38 
14 


-0.4103(6) 

-1.134(7) 

-0.333(1) 


-0.444(1) [21] 

-1.225(2) [21] 



TABLE II. Correlation energies from CIMC with the CCD(l) 
importance function for different r s and N and Ms — 2378. 
For comparison we have also included the results from Ref . [ly 
(basis set independent) and Ref. [2l| (extrapolated using sin- 
gle point extrapolation from M s = 1850 for N = 14 and from 
Ms = 922 for iV = 38, 54). The numbers in parenthesis indi- 
cate the statistical error in each case. 



Slater- Jastrow wave function. These are strict energy 
upper bounds with a bias due to the fixed-node approx- 
imation. Nevertheless, they are believed to be highly 
accurate. 

On the other hand, in Ref. |2l| the energies are calcu- 
lated in a finite CI like basis set (as in this work), us- 
ing the initiator full configuration interaction quantum 
Monte Carlo method (i-FCIQMC). The M s -> oo results 
are obtained by using the so called 'single point extrap- 
olation' from much smaller values of 7V S than ours. 

Our finite-basis set results are already in good agree- 
ment with the other calculations, capturing between 
93% to 99% of the correlation energy. The energy up- 
per bounds can be systematically improved by including 
higher order excitations (triples) in the CC wave function 



and by using larger basis sizes. Possibly, a much faster 
way to achieve basis set convergence is to use a finite 
basis renormalized Coulomb interaction [121 ] . 

Our method shares many similarities with the 
FCIQMC method. The FCIQMC method, in principle, 
can provide exact ground state energies for a CI Hamil- 
tonian. However, most calculations are performed using 
the initiator approximation (i-FCIQMC) with adds a bias 
to the energy estimate. The energies in i-FCIQMC are 
not necessarily upper bounds to the true ground state 
energy. Due to the sign problem, the computational re- 
sources required in either FCIQMC or i-FCIQMC scale 
as the size the many body Hilbert space, i.e., they are 
exponential in the system size, N , and the basis size, M s . 



Computational time (cpu hours) 
N Ms CIMC i-FCIQMC 

200 

2500 

2500 

16000 



TABLE III. Computational cost of our method (CIMC) com- 
pared with the i-FCIQMC method [H| for different r s , N and 

Ms- 



Our method, instead, provides strict energy upper 
bounds, the tightness of which can be systematically im- 
proved by improving the quality of the importance func- 
tion. Importantly, the computational cost of our Monte 
Carlo algorithm nominally grows as N 2 (M S —N), for both 
computational time (per Monte Carlo step) and memory 
requirements. Due to this polynomial scaling, we ex- 
pect our method to be applicable to much larger systems 
than those manageable by conventional diagonalization 
methods or by (i-)FCIQMC. We see evidence of this in 



0.5 


14 


1850 


384 


1.0 


14 


1850 


768 


2.0 


14 


1850 


768 


2.0 


38 


922 


4608 



Table Mil where we compare computational time in our 
method with that in i-FCIQMC. The statistical error for 
the two methods are comparable in each case. 

Conclusion. — We have introduced a many body tech- 
nique which combines the power of ground state projec- 
tion Monte Carlo with coupled cluster theory. As a first 
benchmark and application we studied the homogeneous 
electron gas in momentum space for large single particle 
basis sizes. Our results are already in very good agree- 
ment with existing accurate calculations, and they can 
be systematically improved. In principle, the fixed-node 
bias can be removed by performing a further projection 
starting from our CIMC wave functions with the true pro- 
jector V and signed configurations using FCIQMC and 
the semistochastic projector Monte Carlo method [29]. 

The use of CC wave functions as importance func- 
tions in our Monte Carlo algorithm serves dual pur- 
poses. On the one hand, given that the CC wave func- 
tions are known to be extremely accurate, they can serve 
as the prototype for accurate importance functions in 
Fock space Monte Carlo for normal Fermi systems, anal- 
ogous to the high quality importance functions in r-space 
Monte Carlo |2fl|. 

On the other hand, the energies obtained from conven- 
tional CC theory are not variational. However, the ener- 
gies we obtain from CIMC are rigorous upper bounds for 
the exact ground state energy. By using CC wave func- 
tions in our Monte Carlo we are providing the 'best' vari- 
ational energy one can hope to get using the CC class of 
wave functions. Given the importance of CC type wave 
functions in modern quantum many body calculations, 
the second purpose is as important as the first. 

Computations were performed at the Open Facilities 
at Lawrence Livermore National Laboratory. We would 
like to thank C. Umrigar for stimulating discussions and 
careful reading of the manuscript, and Y. Alhassid for 
comments on the manuscript. 
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